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Abstract 

A new notion is introduced of matrix order indices which relate the matrix 
. norm and its trace. These indices can be defined for any given matrix. They are 

£^ \ especially important for matrices describing many-body systems, equilibrium as 

well as nonequilibrium, for which the indices present a quantitative measure of the 
■ level of ordering. They characterize not only the long-range order, but also mid- 

range order. In the latter case, when order parameters do not exist, the matrix 
indices are well defined, providing an explicit classification of various mid-range 
orders. The matrix order indices are suitable for describing phase transitions 
with both off-diagonal and diagonal order. Contrary to order parameters whose 
correct definition requires the thermodynamic limit, the matrix indices do not 
necessarily need the latter. Because of this, such indices can distinguish between 
different phases of finite systems, thus, allowing for the classification of crossover 
phase transitions. 
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1 Introduction 



Matrices are common objects in quantum physics, where all observable quantities are 
related with matrix representations of operators. In statistical mechanics, a special 
role is played by reduced density matrices [1-3]. The eigenvalues of the latter are 
known to be associated with a possible existence of long-range order called off-diagonal 
long-range order [4-6]. This linkage between the eigenvalues and the order has been 
mathematically formalized by means of order indices for reduced density matrices [7- 
9] . These indices turned out to be a convenient tool for characterizing the arising order 
related to Bose-Einstein condensation and superconductivity [3]. There are, however, 
several limitations not permitting one to employ these indices for describing other 
types of order. Among the most important problems, the following are to be solved: 
(i) How the appearance of magnetic order could be noticed in the behaviour of the 
indices? (ii) How to describe phase transitions associated with diagonal order, such as 
crystallization? (iii) How to characterize crossover-type phase transitions occurring in 
finite systems? 

In the present paper a new notion of matrix order indices is introduced. This notion 
is an essential generalization of the order indices for reduced density matrices [7-9]. 
Contrary to the latter, the matrix indices are defined for any given matrix of arbitrary 
nature. Because of this generality, the matrix indices are suitable for characterizing 
any types of phase transitions, accompanied by any kinds of order, diagonal as well 
as off-diagonal, long-range as well as mid-range, in infinite as well as in finite systems. 
In section 2, a general mathematical definition of the matrix order indices is given. 
These are the quantities linking the norm and the trace of a given matrix whose nature 
can be arbitrary. The possibility of using different norms is discussed. In section 
3, the matrix indices of self-adjoint semipositive matrices are considered, because of 
the special role played by such matrices in physical applications. The convenience of 
employing the introduced matrix indices for quantitatively comparing various kinds of 
mid-range orders is illustrated in section 4. In section 5, Bose-Einstein condensation 
in confined systems is analysed, which is important for describing the condensation 
of trapped atoms. Section 6 demonstrates how the matrix indices must be defined 
in order to treat magnetic phase transitions. Finally, in section 7, it is shown how 
the phase transitions, such as crystallization, characterized by diagonal order, can be 
described by an appropriate definition of the matrix order indices. 



2 Matrix order indices 



Let us have a matrix M of arbitrary nature, with a norm ||M]| and trace TrM. The 
matrix order index is defined as 




(1) 
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We shall say that among two given matrices that one is better ordered whose order 
index is larger. As follows from definition (1), for norm-one and traceless matrices, one 
has ujm = 0. 

The value of the matrix index (1) depends on the type of the norm ||M]| used. 
Assume that a matrix M acts on a normed space {<p} of functions having a vector 
norm \tp\. The matrix norm, associated with a vector norm \tp\, is given by 

„-„ \M<p\ 
M = sup . 

Here and in what follows, it is implied that \ip\ ^ 0. If {ip} is a Hilbert space, with a 
scalar product (<p,(p), the latter naturally generates the Hermitian vector norm \<p\ = 
\J (f, ip). Then, for a matrix M, the Hermitian norm is 

v m 

where M + is a Hermitian conjugated matrix. Supposing that the matrix M possesses 
a spectrm {/x n } of eigenvalues, one has the spectral norm 

\\M\ \ = sup \fj, n \ . (3) 

n 

More generally, the spectral norm is defined as the square root of the largest eigenvalue 
of the matrix MM. When M = M + is a self- adjoint and semipositive matrix, then 
its eigenvalues are semipositive, fj, n > 0. In that case, the spectral norm becomes 

l|M|| = sup^4^ = su P/ , n . (4) 

Among other matrix norms, for physical applications one may employ the trace 
norm 

\\M\\ = (TrM+M) 1/2 , (5) 

which is also called the Euclidean norm. Definition (5) may be convenient for calcu- 
lations, since the trace does not depend on a matrix representation. With a spectrum 
{jjb n } of M, the trace norm (5) reads 



|Mir=(EH 2 ) 



1/2 



As follows from the above definitions, ||M|| < hence one has to distinguish 

between the associated indices % and u' M , for which % < oj' m . 
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3 Self-adjoint semipositive matrices 



Self-adjoint matrices play a special role in physical applications, because of which 
they need to be studied more attentively. Let X = {x} be an arithmetic space of 
physical coordinates and let t G [0, oo) denote time. The set x may include Cartesian 
coordinates, spin, isospin, and other quantum variables. Assume that an algebra of 
local observables A = {A(x, t)} and an algebra of field operators ^ = {ip(x, t),^(x, t)} 
are given on a Fock space T . The direct sum A = A © \P denotes an extended local 
algebra. 

Consider a system of p objects that, for concreteness, will be called particles. The 
set x p = {xi,x 2 , • • • , x p } of particle coordinates pertains to the space X p = X 1 x X 2 x 
. . . x X p , which is a p-fold tenzor product. In a Hilbert space of complex functions 
ip p (x p ), the scalar product is given by the integral 



(<p piV / p ) = J V ;(x p ) <p' p (x p )dx p , 

where dx p = dx\dx 2 ■ ■ ■ dx p . In general, p — 1, 2, . . .. 
Introduce the function 

D p (A,x p ,x p ,t) = TrrA(xi) . . . A(x p )p(t)A + (x p ) . . . , (6) 

where the trace is taken over the Fock space T, A(x) ee A(x, 0) is a member of the 
extended local algebra A, and p(t) is a statistical operator. It is admissible to treat 
the function (6) as an element of the matrix 

D p (A,t) = [D p (A,x p ,x p ,t)} (7) 

with respect to the variables x p . By construction (6), the matrix (7) is self-adjoint and 
positive semidefinite. The action of the matrix (7) on a column ip p ee [(/? p (a; p )] gives the 
column 

D P (A, t)(p p = J D P (A, x p , x p , t) <pp(x p ) dx p . 

A linear envelope of the vectors (p p forms a Hilbert space H p = C{(p p } with a scalar 
product (fpip'p ee ((p p ,(p'p). Each matrix (7) acts on 7i p . The spectral norm for a 
self-adjoint semipositive matrix can be written as 

II n (A + mi ^Dp(A,t)ip p 

\\D P (A, Oil = sup . (8) 

The trace over the variables x p is 

Tr D P (A, t) ee J D P (A, x p , x p , t) dx p . (9) 

According to Eq. (1), one may introduce the matrix indices 

u) P (A,t) ee ^imm , (10 ) 

py ln|Tr D p (A,t)\ 
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For a semipositive matrix, one has 

\\D p (A,t)\\ < Tr D p (A,t) , 

because of which 

u p (A,t)<l. (11) 

Similarly, employing the trace norm (5), one obtains the matrix indices u>'(A,t) for 
which inequality (11) is also valid. 

For constructing a matrix (7), one may choose any representative of the extended 
local algebra A. In particular, opting for the field operators ip(x,t), we come to the 
reduced density matrices 

p p (t) = [p p (x p ,x p ,t)] (12) 

with the elements 

p p {x p ,x p ,t) = D p (ij,x p ,x p ,t) = Trrij(x 1 )...ij(x p ) P (t)ij\x p )...^(x 1 ) . (13) 
The spectral norm for the p-particle density matrix (12) is 

| \p p (t)\ | = sup p pn (t) , (14) 

n 

with {p pn {t)} being the related spectrum. 

The number of particles, N, in a statistical system is assumed to be large. Then, 
in the trace 

AH 

Tr p P {t) = ^— ^ (p= 1,2,. ..,7V) 



one may use the Stirling formula iV! ~ \j2nN N N e N , which, under p fixed and N ^> 1, 
yields 

N\ N p 
(N-p)\ ~ eP 
Hence, the matrix index (10) becomes 

for N > 1. 

For equilibrium statistical systems, one often considers the thermodynamic limit 

N 

iV^OO, V^OO, y^P^ ( 16 ) 

in which is the number of particles in the volume V and the limiting density p is 
finite, < p < oo. In equilibrium, the matrix index (15) does not depend on time, 
becoming 
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as N ^> 1. Note that for calculating the index (17) it is not compulsory to take the 
actual thermodynamic limit (16) by setting N — > oo, but Eq. (17) can be used for 
finite systems, provided that N >> 1. 

From the properties of the reduced density matrices [3], we know that for bosons 
1 1 /5p 1 1 < N p ■ const, hence oj p (ip) < 1, in agreement with the general inequality (11). 
The equality 

u> p (if)) = 1 (boson long — range order) (18) 

corresponds to long-range order in a bosonic system, as that occurring at Bose-Einstein 
condensation of dilute gases. 

Reduced density matrices for fermionic systems [3] satisfy the inequalities 

I \p2pW < N p ■ const , | |p2p-i| | < N p_1 ■ const , 

because of which for the matrix indices we have 

1 p-1 

The limiting case 

^2p(V ; ) — _ (even long — range order) (20) 

happens when even long-range order develops in a fermionic system, like that arising 
in superconductors. When Eq. (20) is valid, then 

" 2hW = ^T (p=l,2,...,JV). 
For large p, this tends to the value (20), since 

uj 2p -iW ~ \ ~ ^ (P> !) ■ 

Even long-range order may also happen for bosonic systems. 

The matrix indices, introduced above, are different from the indices for reduced 
density matrices, defined in the earlier papers [7-9]. The main difference is that, in 
the case of long-range order, the matrix indices oj p (ip) or u> 2p (ip) do not depend on p, 
displaying a kind of resonance, as follows from Eqs. (18) and (20). 

4 Finite momentum condensation 

In order to explicitly illustrate the calculation of the matrix order indices, for both 
long-range and mid-range orders, let us consider a Bose system with finite momentum 
condensate [10-12]. Such a condensate consists of particles having a finite modulus of 
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momentum fc = |k |, with random directions of the vector k . The related condensate 
density matrix reads 

^• r ' ) ^7ik)/ e ' k ' K "' ,rfn(k ° ) ' (21) 

where po = N /V is the condensate density and fi(k ) is a spherical angle associated 
with k . The single-particle density matrix tends, as the distance between the points 
r and r' increases, to the condensate matrix (21), 

Pl (r,r')~p (r,r') (|r - r'| - oo) . (22) 

In a particuar case of zero momentum k = 0, Eq. (21) reduces to the constant density 

Po(r,r')=p (k = 0) 

corresponding to the usual Bose-Einstein condensate with the order index (18). But 
for the finite-momentum condensate, the situation is more elaborate and depends on 
the space dimensionality. 



4.1 Three-dimensional case 

In the three-dimensional space, the condensate matrix (21) reads 

p (r,r') = po 



sin k 


r — r'| 


k 


r 


- r' 





(23) 



For the considered uniform system, the natural single-particle states, which are the 
eigenfunctions of the single-particle density matrix, are the plane waves 



1 



ikr 



Introducing the notation for the radius of the system, 



R = 



3N 



s 47rp) 



1/3 



for the matrix elements of p = [po(r, r')] we have 



(<Pfe,P <Pfe) = 2tt 



PqR 
kk 



sin(/c — k )R sin(/c + k )R 



(k - k )R 



(k + k )R 



In particular, if A: — ^ 0, then 



v / \ a PoR (sin kR N 

hm ( Vk , p oVk ) = 4n — \ - cos kR 



(24) 



From here, the largest eigenvalue of the single-particle density matrix is 



sup lim ((p k ,po<Pk) = N , 

fco— >0 



which results in the order index (18). 
However, for k ^ 0, we find 

sup {(p k ,Po^k) = 2vr — T 

k 



1 - 



sm(2k R) 
2ko~R 



For large systems, such that k R 3> 1, we obtain 



Pi 



2 » % [P\ 



1/3 



which requires that the number of particles be large, 

4np 



3k $ 



Thus, for the order index of the single-particle matrix, we obtain uji(ip) = 1/3. In the 
general case, exploiting the relation 



\Pp\ 



\Pi\ 



valid for asymptotically large N, we have the order indices (17) as 



1 

3 



(25) 



This value of the order indices corresponds to mid-range order. 

If, instead of the spectral norm (4), one uses the trace norm (5), then the order 
indices for reduced density matrices are different. To calculate such indices, we notice 
that 



Trpl = 2W P -f 



o 



sm(2k R) 
2k n R 



For the standard condensate with ko = 0, we have TrpQ = Nq, from where ou[(x/j) = 1, 
which coincides with oJ\{jp) = 1. But if k ^ and the system is large, such that 
koR 1, then we find 

^8 p f3N\ 2/S 



\Pi\ 



k \4irp) 



Finally, we obtain 



(26) 



which differs from the index (25). 
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4.2 Two-dimensional case 

The condensate density matrix (21) is 

Po(r,r / )=p ^o(A:o|r-r / |), (27) 

where 

(-l) n rir . 

J n {z) = ^ \ e tzcoslp cos(n<p) d<p 

IX Jo 

is the first-type Bessel function, with n — 0, 1, 2, . . . Using the relation 

J n (ze™ m ) = e~J n (z) , 
in which m — 0, ±1, ±2, . . ., and the integral 

/x 
xJ n (ax)J n (bx) dx = — — [aJ n+1 (ax)J n (bx) - bJ n (ax)J n+1 (bx)} , 
er — cr 

we find 

(¥k,Po¥k) = ?2 P °l-2 [kJ\{kR)Jo{k Q R) - k Q J Q (kR)Ji(k Q R)} , 

where 

'N\ 1/2 



is the system radius. 
When k — > 0, then 



Rs [ir P (28) 



PqR 

]im (<p k , po<p k ) = 2tt -5— J^kR) . 

k ^0 K 



Because of the asymptotic equation 

1 (z\ n 



it follows that 



sup27r ^— JAkR) = 7iR 2 p = N 
k k 



And we again return to the indices (18). 

However, for the finite-momentum condensate, with k ^ 0, we have 



SUp ((fik:P0<Pk) = np R 2 

k 



Jl{k Q R) - MhQMhR) + jKkoR) 



where the equations 

dJ n (z) 
dz 



nJ n (z) - zJ n+ i(z) = zJ n -i(z) - nJ n (', 
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have been used. For a large system, such that koR 3> 1, we invoke the asymptotic form 

J n(z) - \[^ cos (z - | n - ^ (z->oo). 
Then, for a system of many particles, with 



Tip 

iV 

we get 



||pi||-4^f^y /2 cos 2 (A; i?). 
For the order indices we obtain 

= \ . (29) 

If calculating the order indices, one opts for the trace norm (5), then the resulting 
u' p (ip) differs from the indices (29). Thus, using the relations 

j. n (z) = (-iruz) , 

X.J^(ax) dx = Jn( aX ) ~ Jn-l{ ax )Jn+l{Q>x) 

2 

we have 

Tr pi = ( Po nR 2 ) [j%{k R) + J*{k R) . 

The indices oo' (ip) coincide with uj p (ip) = 1 only if /c = 0, since then Tr/?Q = N$. 
But for k ^ 0, again keeping in mind a large system, with k$R ^> 1, we find 

[T/N\ V4 

m ^ 2p ik\T P ) i cos (^)i. 

As a result, we obtain 

= \ , (30) 

which is different from the indices (29). 

4.3 One-dimensional case 

The condensate density matrix (21) takes the form 

p (x,x') = p cosk (x - x') , (31) 

where x G (— oo, +oo). With the system length 

N 

L = - , (32) 
P 
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we have 



sin L) sin (^a ^ 



+ 



(A; - A; )L (A; + fc )L 
If fc = 0, then again we get oo p (ip) = 1, since 

lim o (^,po^) = 2f sin 

and the supremum over k for the above right-hand side gives N Q . 
The situation is rather interesting for k ^ 0. Then 



SUp {(pk,P0<Pk) = -7T- 



1 + 



sin(A;oL) 
k L 



For a large system, with k Q L ^> 1, and, respectively, with a large number of particles, 
such that 



we find 



Pi 



Po 
2p 



AT . 



Finally, we obtain 

u p (if>) = 1 . 

When employing the trace norm (5), and using the equality 

sin koL \ 



(33) 



^Po 



\(PoL) 2 (l 



k L 



we get for a large system 
This yields 



i . 



(34) 



Hence, in the case of one-dimensional systems, the order indices (33), based on the 
spectral norm, coincide with the indices (34), based on the trace norm of reduced 
density matrices. 

The order indices, defined through the spectral norm, or respectively through the 
trace norm, permit us to generalize the consideration to an arbitrary <i-dimensional 
case resulting in 

= \ , = ^ ■ (35) 



2d 
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This suggests a relation 

between the order indices obtained for different norms. 

Let us emphasize that the order indices provide us with a quantitative measure for 
the degree of order occurring in a statistical system. This measure makes it possible 
to unambiguously distinguish between the degrees of order of quite similar systems. 
For example, let us compare two three-dimensional systems, one with the condensate 
density matrix (23) and another with that matrix having the same algebraic decay, but 
without oscillations due to sine, that is when p (r,r') ~ p /ko\r — r'|. For the latter 
case, we find u> p (ip) = 2/3 instead of the index (25), which tells us that in the latter 
system the order is higher than in the system with the condensate matrix (23). 

It is important that the order indices allow the description and classification of 
various types of mid-range order, as is illustrated in this section, when the order pa- 
rameters are zero. Note that a Bose system with finite- momentum condensate [10-12], 
possessing mid-range order, satisfies the Cummings-Hyland- Rowlands relation [13]. 

5 Condensation in confined systems 

The recent developments of Bose-Einstein condensation of trapped atoms (see reviews 
[14-16]) makes it important to consider how the order indices can be defined for finite 
systems. Stricktly speaking, in finite systems sharp phase transitions cannot occur. 
A rigorous definition of a genuine phase transition requires the usage of the thermo- 
dynamic limit (16). But the number of atoms in a trap, although large, however is 
always finite. Therefore the Bose-Einstein condensation in a trap is rather a gradual 
crossover than a sharp transition. The introduced matrix order indices do not neces- 
sarily require the employment of the thermodynamic limit, though their calculation is 
facilitated when the number of particles is large. 

Bose-Einstein condensation can be treated as the appearance of a macroscopic 
fraction of atoms being in a coherent state. A coherent state, in general, is defined up 
to a phase factor with a random phase [17]. If this phase is made fixed, then the gauge 
symmetry in the system is broken. The breaking of gauge symmetry is usually done by 
means of the Bogolubov prescription [18] presenting the field operator as ip — i/jq + i/j, 
where ip is a classical quantity and ip is the field operator of noncondensed particles. In 
that case, the statistical average < ip >= ip is nonzero and might be considered as an 
order parameter. The gauge symmetry breaking implies that the number of particles 
in the system is not conserved. In reality, the number of particles in a system, even 
such as a trap, can be preserved with a rather good accuracy. Consequently, it would 
be logical to describe Bose-Einstein condensation without breaking gauge symmetry 
[19-24]. The suggested schemes, preserving particle conservation, are essentially more 
cumbersome then the techniques based on gauge symmetry breaking. This is why 
the latter techniques became more popular. However this way also contains intrinsic 
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difficulties related with the appearance of divergences in low-order approximations. 
To compensate these divergences, one needs to invoke higher-order approximations 
[23], which again complicates the consideration. Here we advance a novel method 
of treating a system with Bose condensate, without breaking gauge symmetry. The 
method is sufficiently simple and general, and can be applied to trapped as well as 
untrapped atoms interacting through integrable potentials. 

The pivotal idea of the method is the usage of random-phase coherent states [17]. 
When a fraction of atoms is in such a state, the field operator can be presented as the 
sum 

iP(r,t) =r} a (r,t)+$(r,t) , (36) 

in which 

r ta (r,t)=r t (r,t)e ia (37) 

is a random-phase coherent field [17] with a G [0, 2n] being a random phase, and ip 
corresponds to incoherent particles. 

The average of an operator A is defined as 

1 r 2n 

< A > = — / TrpA da , (38) 



2tt Jo 

including statistical averaging, with a statistical operator p, and averaging over the 
random phases. From this definition, it follows that 

<^> = <^ a > = <^> = 0, (39) 

hence the gauge symmetry is not broken. The conservation of gauge symmetry is 
what makes the principal distinction between the presentation (36) and the Bogolubov 
prescription [18], where the gauge symmetry is broken. The coherent and incoherent 
terms in the sum (36) are subject to the normalization conditions 

J <^(r,t)i>(r,t)> dr = N, J \r)(r, t)\ 2 dr = N , 

< ft(r,t)i>(r,t) > dr = N , (40) 



/ 



where N is the number of coherent atoms, N is the number of incoherent atoms, and 
the total number of atoms is 

N = N + N . (41) 

The terms in the right-hand side of Eq. (36) are mutually orthogonal on the average, 
which means that 

<r] a (r,t)4>(r',t) > = 0. 

For the standard Hamiltonian with a two-particle interaction potential $(r) = 
$(— r), the Heisenberg equation for the field operator reads 

ih ^ = , (42) 



13 



where tp = tp(r,t) and 



Hty) = ~ ^ + U ext + U mt - p , (43) 

with m being particle mass and p, chemical potential. The external potential U ext (j, t) 
describes all external fields, including the confining potential and possible external 
perturbations. The interaction potential 

U int (r, t) = J $(r - r')^(r', t)i>(r', t) dv' (44) 

is formed by particle interactions. 

Let us substitute the presentation (36) into the Heisenberg equation (42). Then, 
let us either multiply Eq. (42) by e~ ia or just set a = 0, after which take the average 
< . . . > according to definition (38). This results in the equation for the coherent field, 

(^ + — -U ext + p) V (r,t) = 

= J $(r-r') [\r ) (r',t)\ 2 r ) (r,t)+p(r',t)r ] (r,t) + p 1 (r,r',t)r)(r',t)\ dv' , (45) 
in which the density of incoherent particles 

p(r,t) =pi(r,r,f) (46) 
is the diagonal element of the incoherent density matrix 

p 1 (r,r',t) = <^(r',t)^(r,t)> . (47) 

Notice that in the equation (45) for the coherent field no anomalous averages arise, 
which is due to the gauge symmetry conservation. 

If, after substituting the presentation (36) into the Heisenberg equation (42), we 
average only over the random phase, then we obtain the following equation 

/ d U 2 V 2 \ ~ 

(^ + — -U ext + ^^r,t) = 

= J ^(r - O [^(r', ^^(r/, «)^( r> + ^(r', «)|^( r , t) + » 7 *(r', «)»7(r, *)^(^» *)] dr> 

(48) 

for the incoherent term of the field operator. Let us stress that both equations (45) as 
well as (48) are exact. 

One often obtains an equation for the condensate wave function, which is equivalent 
to Eq. (45), following what one calls the Popov approximation. In that way, one, first, 
breaks the gauge symmetry by means of the Bogolubov prescription and then one omits 
the anomalous averages arising because of the broken gauge symmetry. Such a way is 
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clearly not self-consistent. Moreover, as is possible to show by direct calculations, in 
the case of broken gauge symmetry, anomalous averages are of the same order, or even 
larger, then normal averages [25,26]. Therefore, after the gauge symmetry has been 
broken, it is not correct to omit anomalous averages. 

Given Eq. (48), we may derive an equation for the Green functions of incoherent 
particles. For this purpose, it is convenient to introduce the following short-hand 
notation [17]. Denote a function /(r^ti) as /(l) and the related differential measure 
as d(l) = dr-idti. Then the Dirac delta-function is 

5(12) = <S(n - r 2 )<S(*i - h) . 

Introduce the retarded interaction potential 

$(12) = $(ri-r 2 )<J(*i-* 2 + 0) . 

With this notation, the single-particle and two-particle causal Green functions are 
defined as 

G(12) = -i < f ^(1)^(2) > , G(1234) = - < f ^(1)^(2)^(3)^(4) > , (49) 

where T is the chronological operator. These are the Green functions at real times, 
which are directly related to the Matsubara Green functions at imaginary times [17,27]. 
From definition (49) and Eq. (48), we immediately obtain the propagator equation 

- J $(13) [|r ? (3)| 2 G(12) + 77*(3)?7(1)G(32) + iG 2 (1332)] d(3) = 5(12) . (50) 

It is important to note that all Eqs. (45), (48), and (50) have sense only for the 
integrable interaction potentials, such that 



$(r) dr 



< oo . (51) 



In the other case, to avoid ultraviolet divergences, one has to set i](r,t) = 0, which 
implies the absence of condensate [17]. 

The propagator equation (50) can be written in a more compact form if one intro- 
duces the self-energy £(12) by means of the relation 



£(13)G(23) d(3) = 

= J $(13) [\r](3)\ 2 G(12) + rj*(3)rj(l)G(32) + iG 2 (1332)] d(3) . (52) 
From here, recalling the definition of the inverse propagator, 

J G- 1 (13)G'(32) d(3) = J G(13)G- 1 (32) d(3) = 5(12) , 
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we have the self-energy 

S(12)=5(12)|$(13)|^(3)| 2 rf(3) + 
$(12)?7(1)77*(2) +i J $(13)G 2 (1334)G'- 1 (42) d(34) . 



+ 



Then the inverse propagator is 



^(12) = (in^ + g-^ + ^(12)- E (12) 



(53) 



(54) 



and the propagator equation (50) takes the form 
d Ti 2 V 2 



ih 



^ + ^ - U ext + t?j G(12) - J S(13)G(32) d(3) = 5(12) 



(55) 



The pair of equations (45) and (55), together with the definition (53) for the incoherent 
self-energy, completely describes the Bose system with a coherent condensate. These 
equations, for a given interaction potential, are exact. 

Dilute gases of trapped atoms are usually characterized by the Fermi contact po- 
tential 



$(r) = AS(r) , 



A = 4nh 2 ^ , 
m 



(56) 



where a s is a scattering length. 

In that case, the coherent-field equation (45) simplifies to 



ih 



dr] 
dt 



h 2 v 2 

2m 



+ U ext + A (\r]\ 2 + 2p) - /I 



(57) 



where rj = r](r,t) and p = p(r,t). Equation (48) for the incoherent operator becomes 

h 2 v 2 



dip 

th m = 



2m 



+ U ext + A + 2\r]\ 2 ) - n 



The self-energy (53) transforms to 



£(12) = A 



25(12)|??(l)| 2 + i J G 2 (1113)G'- 1 (32) d(3) 



(58) 



(59) 



which enters the propagator equation (55). 

To solve the propagator equation, one may follow one of the known iterative pro- 
cedures [17,28] starting e.g. with the Hartree-Fock approximation, when 

G 2 (1234) = G(14)G(23) + G(13)G(24) . 

Then the self-energy (59) is 

£(12) = 2A5(12) [|??(l)| 2 + p(l)l . (60) 
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If the external potential does not depend on time, U ext = U(r), and the system is 
in equilibrium, then the coherent field may be presented as 

*7(r, t) = y/N~ <p(r) exp (- % - Etj , (61) 

with the function <p(r) normalized to unity, (<p, ip) = 1. In this case, Eq. (57) reduces 
to the stationary equation 

H(p)p(v) = Ep(v) , (62) 

with the effective Hamiltonian 

H(p) = - T ^ + U + A (n \^\ 2 + 2p) -ii. (63) 

Note that this Hamiltonian depends on thermodynamic parameters, such as tempera- 
ture, through the incoherent density p = p(r). 

In the presence of a confining potential, the spectrum of the eigenproblem (62) is 
discrete [29]. For a system in equilibrium, the condensate energy E corresponds to 
the lowest energy level of the eigenproblem (62). For a nonequilibrium system, one 
can excite higher nonlinear coherent modes corresponding to higher states of Eq. (62), 
thus, creating nonground-state coherent condensates [29-33]. 

For an equilibrium system, the single-particle density matrix 



is the sum 

of the coherent term 



Pi(r, r') = < ip'(r')ip(r) > 
p 1 (r,r')=p (r,r / )+Pi(r,r / ) (64) 



p (r,r') = N p (r)Po( r ') , 

where <po( r ) is the ground state of the eigenproblem (62), and of the incoherent density 
matrix (47). The coherent density matrix p = [po(r, r')] has the properties 

p^ = <- 1 po, Trpg = <. 
The incoherent density matrix can be presented as an expansion 

Pi(r,r') =J2n k Pk(r)pl(r') 

k 

over natural wave functions [3] corresponding to incoherent states. 

When the coherent condensate is absent, iV = 0, then the spectral norm of pi is 



||pi|| = sup k nk, which is finite. The trace norm gives ||pi||' = yJ2k n li which is of 
order N 1 ^ 2 . The related order indices are uj p (vp) = and uo' p (ip) = 1/2, respectively. 
However, as soon as iVo > sup fc nfc, then for a large system with N ^> 1, one has 

MVO = = ^ • (65) 
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In the thermodynamic limit (16), we return to the order indices (18). But for a finite 
system the order indices (65) are, in general, functions of thermodynamic variables, 
such as temperature and density. When N changes from 1 to N, the order indices 
increase from to 1, which corresponds to a gradual crossover. 

Recall that the present consideration has been based on the decomposition (36) 
for the field operator onto its coherent and incoherent parts. This decomposition is 
valid only when the interaction potential is integrable, satisfying condition (51). The 
Fermi contact potential (56) is integrable. This potential describes pair interactions 
in dilute gases, whose mean interatomic distance is much larger than the magnitude 
of the scattering length \a s \. For dense gases and liquids, interaction potentials often 
contain hard cores and are not integrable. In such cases, the decomposition (36) is 
not appropriate and one has to invoke other calculational techniques explicitly taking 
account of interatomic correlations (see e.g. [34-36] and references in review [16]). 



6 Magnetic phase transitions 

The matrix indices (17), defined for reduced density matrices, are convenient for char- 
acterizing such phase transitions as Bose condensation and superconductivity. The 
question that one could pose is whether the same indices (17) would be suitable for 
describing magnetic phase transitions. And also, since superconductivity and magnetic 
order can coexist [37,39], could the indices (17) characterize such a coexistence of or- 
ders? In this section, we show that the order indices (17), based on reduced density 
matrices (18), are not suitable for describing magnetic order. To describe the latter, it 
is necessary to return to the general definition (1) of the matrix order indices and to 
introduce the matrices composed of spin operators. 

Consider, first, the indices u} p (ip). Let the field operator be a spinor ip = with 
a =|, | denoting spin up or down. Keeping in mind a solid with a crystalline lattice, 
we enumerate the lattice sites by the index % = 1, 2, . . . , N. Wave functions, localized 
at the related sites, are called the localized orbitals, such as Wannier functions. For 
simplicity, we shall deal with the single- zone case. Let {</?j(r)} be a set of localized 
orbitals that are orthonormalized, (ipi,ipj) = 5ij, and asymptotically complete, which 
means that 

N 

E^(rM(r')=J N (r-r'), (66) 
with the right-hand side asymptotically, as iV — > oo, coinciding with the 5-function, 

lim Mr) = *(r) ■ (67) 
The field operator ip a (r) can be expanded over the localized orbitals as 

N 

Vv(r) = ]>>^(r) . (68) 
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Both operators ip a ( r ) as well as q ct satisfy the Fermi commutation relations. To stress 
that the particles are well localized, the following localization conditions are imposed 
on the operators Cj CT : 

Limiting the number of particles at each site by one, we have the unipolarity conditions 

J2 C L C i« = 1 , °iaCia' = . (70) 

From the particle operators one can pass to spin operators S" by means of the 
transformations s _ 

c lt c *T = 2 ~*~ , c li c il — 2 — ' 

c 1t c *I = 3? = + > c ij. c «T = = &i ~ • 

Now, we construct the reduced density matrices (13). The single-particle density 
matrix is 

Pl (r,r') = Mr-r'), (71) 

where Sn{^) is the asymptotic 5-function, appearing in Eq. (66). Note that, because 
of Trpi = N, one has 5n(0) = N/V . The spectral and trace norms, respectively, give 

||pi||~l, ||pi||' ~ VN (iV->oo). 

From here 

WiM=0, wlW = \, (72) 

independently from the existence or absence of magnetic order. 
For the two-particle reduced density matrix, we have 

p 2 (ri, r 2 , r;, r' 2 ) = ^(n, ri)pi(r 2 , r' 2 ) - ^ ^(n, r^Pifo, r'J- 

iV 

- 2]T < St ■ Sj > ^(rO^-Cra)^^)^^) . (73) 

To simplify the consideration, let the particle interactions be of long-range type, for 
which the asymptotic, as iV — > oo, decoupling 

< Si • Sj > < Si >< Sj > (i^ j) 

is valid. And also let the lattice be ideal, when, because of the translation invariance, 
one has 

1 N 

<s t > = <s>, s = -J2s t . 

i=i 
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Then the two-particle matrix (73) becomes 

p 2 (n, r 2 , ri, r 2 ) = Pl (n, ri)pi(r 2 , r' 2 ) - Q + 2 < S > 2 ) ^(n, r 2 ) Pl (r 2 , ri) . (74) 
The eigenfunctions of this matrix are 

fij( r , r ') = ^= [^(r)^(r') - ^(r')^(r)] , 
with (fii(r) being the localized orbitals. For a large system, we get 

wfow ~ ^+2 < s > 2 , Hp2ir~iv, 

as iV — > oo. Thus, in the thermodynamic limit, we find 

u> p ty) = , u/ p ty) = i • (75) 

This shows that, although the norms of the reduced density matrices do change un- 
der the arising magnetic order, that change is negligible in the thermodynamic limit. 
Hence, the indices uj p (ip) are not suitable for describing magnetic phase transitions. 

This problem can easily be resolved by going back to the more general definition 
of the matrix indices (10). Being interested in magnetic order, it is reasonable to 
work with correlation functions composed of spin operators. For instance, the two-spin 
correlator can be written as 

D ij (S) = <S j 'S i > . (76) 

In general, the 2p-spin correlator is 

D h .., pjl ... jp (S) =< K ■ ■ ■ (S n • S n ) • • • S ip > . (77) 

The related spin- correlation matrices D P (S) are the matrices with the elements (77). 
For example, Di(S) = [Dij(S)] is an N x N matrix with the elements (76). 
For a system defined on an ideal lattice, the natural wave functions are 

M*) = 77= ^ , (78) 



with aj being a lattice vector. These functions form an orthonormalized basis, for 
which 

i k 

Calculating the spectral norm of Di(S), we find 

JV-l 



1^(5)11 = sup 

k 



s(s + 1) + J2 < s j ■ s i > e 

3=1 



ika, 
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The trace norm yields 

(N-l 
NJ2< S, • So > 

And for the trace, we have 

Tr Di(S) =S(S+1)N . 
When magnetic order is absent, i.e. < S >= 0, then for a large system, with N 1, 

||A(5)||~1, WhfflW' ~ N 1 ' 2 . 
Hence, generalizing this, we get 

u p (S)=0, <J p {S) = \ (<S> = 0). (79) 

But as soon as there appear magnetic order, then 

\\Di{S)\\ ~N < S > 2 , \\Di{S)\\' ~ N\ < S > | , 

as iV — > oo. In this way, we come to the order indices 

u p (S) = u' p (S) = 1 (< S > ^ 0) . (80) 

Comparing Eqs. (79) and (80), we see that the order indices w p (S) perfectly distinguish 
between the existence and absence of magnetic order. The description of the latter 
by the indices cu p (S) is analogous to the characterization of Bose condensation by the 
indices oj p {ip). This is not surprising, since the arising magnetic order can be interpreted 
as condensation of magnons [39]. 

7 Transitions with diagonal order 

The eigenvalues of the reduced density matrices p p are known to be sensitive to the 
existence of the off-diagonal long-range order [3-6]. This is why the indices oj p {jp) are 
so useful for characterizing phase transitions with arising off-diagonal order, such as 
superconductivity and Bose condensation. Similarly, the indices cu p (S) are convenient 
for describing magnetic transitions. These indices, however, are not sensitive to the 
appearance of the so-called diagonal long-range order, associated with the properties 
of the diagonal elements of density matrices. Among phase transitions with developing 
diagonal order, we may mention the melting-crystallization phase transition, which is 
so well known on earth and, presumably, occurring in neutron stars [40]. The mixing- 
stratification transformation is also of this kind. 

The difference between phase transitions with off-diagonal long-range order and 
diagonal long-range order is in the following. In the former case, the largest eigenvalue 
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of one of the reduced density matrices, and respectively the related spectral norm, 
change from being of order one to the order of N p . Such a dramatic variation of the 
norm, due to the arising off-diagonal order, immediately results in the sharp change 
of the order index w p (ip). But so drastic increase of eigenvalues of the reduced density 
matrices does not happen under the appearing diagonal order, when all eigenvalues of 
p p remain of order one. Thus, for both crystalline and liquid states 

Therefore the indices 

uM) = , u£M = \ (81) 

do not notice that the diagonal order has appeared. 

Fortunately, with the general definition of the matrix order indices (1) or (10), we 
are not obliged to deal exclusively with the density matrices (12). It is always possible 
to define such a matrix whose order index would be sensitive to any kind of arising 
order. 

As an illustration, let us consider the crystallization transition. Introduce the 

density- difference operator 

A(r) = ^(r)^(r) - p , (82) 

with p = N/V being the average density. Statistical averaging of the operator (82) 
gives 

< A(r) > = p(r) - p , p(r) = < ^(r)^(r) > . 
Define the matrix -Di(A) whose elements are 

Di(A, r, r') = < A(r')A(r) > . (83) 

Similarly, we may construct a matrix D p (A) of any order p — 1, 2, . . ., as is explained 
in Eq. (6). The elements of D P (A) satisfy the correlation weakening condition, e.g. 

< A(r)A(r') > ~ < A(r) > < A(r') > (|r - r'| -> oo) . 

In the case of liquid, with uniform density, we have 

< A(r) > = , p(r) = p . 

Eigenvalues of Di(A) are of order 1, which yields 

||A(A)||~1, WD^AW ~ N 1 ' 2 . 

Following this way, we find 

u p (A) = , 4(A) = \ , (84) 
which signifies the absence of order. 
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The occurrence of crystalline order is marked by a well defined crystalline lat- 
tice spanned by the lattice vectors a^, with i = 1, 2, . . . , N. Introduce the function 
D 1 (A,a i ,a.j) as in Eq. (83) and the related matrix -Di(A) = [-Di(A, a*, a,-)], which is 
an N x N matrix. The natural wave functions for a lattice system are presented in Eq. 
(78). Employing the spectral norm, we get 

||A(A)||~7V[p(a)-p] 2 (N^oc), 

where a is any of a i; since for an ideal lattice, p(aj) = p(sLj). The same expression 
follows for the trace norm ||£) 1 (A)|| / . Finally, we obtain 

u p (A) = u' p (A) = 1 . (85) 

As is seen, the formation of crystalline order leads to the change of the matrix indices 
LOp(A) from the values in Eq. (84) to those in Eq. (85). 

As s brief conclusion, let us stress that the matrix order indices, introduced in the 
present paper, is a new general notion allowing for the description and classification 
of arbitrary types of order, off-diagonal as well as diagonal, long-range and mid-range, 
occurring in infinite as well as finite systems, and happening in different spaces (real, 
momentum, spin,. . .). The matrix order indices can be well defined when the order pa- 
rameters do not exist. These indices are suitable for describing sharp phase transitions 
as well as gradual crossovers. 
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